An explicit GPU-based material point method solver for elastoplastic problems (ep2-3De v1.0)

We propose an explicit GPU-based solver within the material point method (MPM) framework using graphics processing units (GPUs) to resolve elastoplastic problems under two- and three-dimensional configurations (i.e. granular collapses and slumping mechanics). Modern GPU architectures, including Ampere, Turing and Volta, provide a computational framework that is well suited to the locality of the material point method in view of high-performance computing. For intense and non-local computational aspects (i.e. the back-and-forth mapping between the nodes of the background mesh and the material points), we use straightforward atomic operations (the scattering paradigm). We select the generalized interpolation material point method (GIMPM) to resolve the cell-crossing error, which typically arises in the original MPM, because of the C0 continuity of the linear basis function. We validate our GPU-based in-house solver by comparing numerical results for granular collapses with the available experimental data sets. Good agreement is found between the numerical results and experimental results for the free surface and failure surface. We further evaluate the performance of our GPU-based implementation for the three-dimensional elastoplastic slumping mechanics problem. We report (i) a maximum 200-fold performance gain between a CPU- and a single-GPU-based implementation, provided that (ii) the hardware limit (i.e. the peak memory bandwidth) of the device is reached. Furthermore, our multi-GPU implementation can resolve models with nearly a billion material points. We finally showcase an application to slumping mechanics and demonstrate the importance of a three-dimensional configuration coupled with heterogeneous properties to resolve complex material behaviour.</p

Resolving Wave Propagation in Anisotropic Poroelastic Media Using Graphical Processing Units (GPUs)

Biot's equations describe the physics of hydromechanically coupled systems establishing the widely recognized theory of poroelasticity. This theory has a broad range of applications in Earth and biological sciences as well as in engineering. The numerical solution of Biot's equations is challenging because wave propagation and fluid pressure diffusion processes occur simultaneously but feature very different characteristic time scales. Analogous to geophysical data acquisition, high resolution and three dimensional numerical experiments lately redefined state of the art. Tackling high spatial and temporal resolution requires a high-performance computing approach. We developed a multi- graphical processing units (GPU) numerical application to resolve the anisotropic elastodynamic Biot's equations that relies on a conservative numerical scheme to simulate, in a few seconds, wave fields for spatial domains involving more than 1.5 billion grid cells. We present a comprehensive dimensional analysis reducing the number of material parameters needed for the numerical experiments from ten to four. Furthermore, the dimensional analysis emphasizes the key material parameters governing the physics of wave propagation in poroelastic media. We perform a dispersion analysis as function of dimensionless parameters leading to simple and transparent dispersion relations. We then benchmark our numerical solution against an analytical plane wave solution. Finally, we present several numerical modeling experiments, including a three-dimensional simulation of fluid injection into a poroelastic medium. We provide the Matlab, symbolic Maple, and GPU CUDA C routines to reproduce the main presented results. The high efficiency of our numerical implementation makes it readily usable to investigate three-dimensional and high-resolution scenarios of practical applications.ISSN:2169-9313ISSN:0148-0227ISSN:2169-935

Resolving coupled physical processes in porous rocks: From linear quasi-static and dynamic phenomena to non-linear instabilities

La majorité des processus physiques sur Terre sont couplés. Un processus physique peut en induire un autre, ce qui est le cas d’une onde se propageant dans une roche poreuse saturée de fluide, ce qui induit un écoulement de fluide. Dans un milieu biphasique, l’interaction entre les phases solides et fluides conduit à des effets physiques qui ne sont pas observés dans un milieu monophasique. Ainsi, une description satisfaisante de systèmes physiques complexes nécessite un traitement particulier. L’application des théories décrivant les processus physiques couplés dans les roches poreuses fracturées est d’une grande importance dans les scénarios impliquant la séquestration géologique du CO2, l’élimination des déchets nucléaires, l’exploration et la production d’énergie géothermique et les hydrocarbures. Le développe- ment de méthodes géophysiques non invasives de détection et de surveillance de ces formations géologiques est crucial. La recherche scientifique vise à faire progresser la description quantitative et qualitative des processus physiques couplés dans les roches poreuses. Conformément à ces objectifs, les contributions présentées ici sont réparties dans quatre disciplines différentes : la micromécanique, la géophysique, la mécanique computationnelle et la poroélasticité computationnelle, et la théorie des instabilités non linéaires. Différentes méthodes analytiques et numériques sont utilisées pour résoudre la physique aux micro- et macro-échelles. Cela comprend l’étude des processus linéaires quasi-statiques et dynamiques. Ce travail de recherche contient également des résultats basés sur la théorie des processus physiques non linéaires. À l’échelle microscopique (ou à l’échelle des pores), une roche est constituée de grains, de pores et de fractures individuels. De petites déformations causées par une onde sismique se propageant à travers la roche induisent des gradients de pression dans les fractures conformes. En conséquence, un écoulement de fluide (appelé “écoulement de jet”) a lieu jusqu’à ce que la pression interstitielle s’équilibre. En raison de la viscosité du fluide et du frottement visqueux associé, un tel écoulement de fluide provoque une forte dissipation d’énergie. Des simulations numériques tridimensionnelles de l’écoulement de jet à l’aide d’une approche par éléments finis sont effectuées et les résultats sont comparés à un modèle analytique publié. A partir de cette comparaison, de nombreuses limitations de la solution analytique publiée sont quantifiées et décrites. Par la suite, un nouveau modèle analytique pour la dispersion sismique et l’atténuation associées à l’écoulement de jet est présenté pour une géométrie de pores qui a été classiquement utilisée pour expliquer ce mécanisme d’expulsion de l’eau. Ensuite, ce modèle analytique est étendu pour traiter des géométries plus complexes de l’espace poreux, beaucoup plus représentatives de celle d’une roche. Les paramètres clés de l’espace poreux qui contrôlent la fréquence caractéristique (à laquelle se produit le maximum d’atténuation) sont redéfinis. De plus, des expressions analytiques pour calculer les propriétés de rigidité effective d’un modèle de roche, dont l’espace poreux est décrit par une fracture reliée à un pore ou à plusieurs pores, sont fournies. À l’échelle macroscopique, une roche poreuse peut être décrite par un ensemble de propriétés macro- scopiques, par exemple, les modules élastiques effectifs, la perméabilité, etc. Les équations de Biot décrivent un système couplé hydro-mécaniquement et établissent la théorie largement reconnue de la poroélasticité. Le milieu biphasique est représenté par une matrice poreuse solide élastique saturée d’un fluide visqueux compressible. La réponse dynamique d’un tel milieu biphasique et isotrope se traduit par deux ondes longitudinales et une onde de cisaillement, comme prédit par Yakov Frenkel. La modélisation numérique efficace et précise des équations de la poroélasticité de Biot nécessite la connaissance des conditions exactes de stabilité. Cette recherche présente les résultats de l’analyse de stabilité de von Neumann des équations de Biot discrétisées et de l’équation d’onde amortie linéaire discrétisée. Les conditions exactes de stabilité pour un certain nombre de schémas implicites et explicites sont dérivées. De plus, un solveur numérique d’unités de traitement multi-graphiques (GPU) est développé pour résoudre les équations élastodynamiques anisotropes de Biot afin de simuler, en quelques secondes, des champs d’ondes pour des domaines spatiaux impliquant plus de 4,5 milliards de mailles. Une analyse dimensionnelle complète est présentée, réduisant ainsi le nombre de paramètres matériels nécessaires pour les expériences numériques de dix à quatre. Une analyse de dispersion en fonction de paramètres adimensionnels est effectuée, conduisant à des relations de dispersion simples et transparentes. La haute efficacité de notre implémentation numérique la rend facilement accessible pour étudier des scénarios tridimensionnels et à haute résolution d’applications pratiques. Dans le cadre de la théorie des instabilités non linéaires, une nouvelle théorie de la nucléation sismique est présentée. La rhéologie visco-plastique ou élasto-plastique la plus simple permet de modéliser la nucléation sismique spontanée. En augmentant lentement la contrainte dans le milieu, elle atteint la limite plastique, produisant ainsi la localisation de la déformation et entraînant le développement lent de bandes de cisaillement fractales. Au fil du temps, ces dernières se développent spontanément et des chutes de contrainte se produisent dans le milieu. Une chute de contrainte correspond à un nouveau modèle particulier de localisation de déformation, qui agit alors comme source sismique et déclenche la propagation des ondes sismiques. Cette nouvelle approche de modélisation est basée sur des lois de conservation sans aucune relation constitutive dérivée expérimentalement. -- The majority of the physical processes on the Earth are coupled. A physical process might induce a different one, which is the case of a wave propagating in a fluid-saturated porous rock and inducing fluid flow. In a two-phase medium, the interaction between solid and fluid phases leads to physical effects, that are not observed in a single-phase medium. Thus, a successful description of complex physical systems requires special treatment. The applications of theories describing coupled physical processes in cracked porous rocks are of great importance in scenarios involving CO2 geological sequestration, nuclear waste disposal, the exploration and production of geothermal energy, and hydrocarbons. Developing non-invasive geophysical detection and monitoring methods for these geological formations is crucial. Scientific research aims to advance the quantitative and qualitative description of coupled physical processes in porous rocks. In line with these objectives, the contributions presented here are distributed across four different disciplines: micromechanics, geophysics, computational mechanics and computational poroelasticity, and the theory of non-linear instabilities. Different analytical and numerical methods are used to resolve the physics at the micro- and macro-scales. It includes the study of linear quasi-static and dynamic processes. This research work also contains some results based on the theory of non-linear physical processes. At the micro-scale (or pore-scale), a rock consists of individual grains, pores, and cracks. Small deformations caused by a passing seismic wave propagating through the rock induce pressure gradients in compliant cracks. As a result, fluid flow (so-called squirt flow) takes place until the pore pressure equilibrates. Due to the fluid viscosity and the associated viscous friction, such fluid flow causes strong energy dissipation. Three- dimensional numerical simulations of squirt flow using a finite-element approach are performed and the results are compared against a published analytical model. From this comparison, many limitations of the published analytical solution are quantified and described. Subsequently, a new analytical model for squirt flow associated seismic dispersion and attenuation is presented for a pore geometry that has been classically used to explain squirt flow. Then, this analytical model is extended to deal with more complex geometries of the pore space, which are much more closely representative of that of a rock. The key parameters of the pore space which control the characteristic frequency (at which the maximum of attenuation occurs) are re-defined. Additionally, closed-form analytical expressions to calculate the effective stiffness properties of a rock model whose pore space is described by a crack connected to a pore or multiple pores are provided. At the macro-scale, a porous rock can be described by a set of macroscopic properties, e.g., effective elastic moduli, permeability, etc. Biot’s equations describe a hydro-mechanically coupled system and establish the widely recognized theory of poroelasticity. The two-phase medium is represented by an elastic solid porous matrix saturated with a compressible viscous fluid. The dynamic response of such an isotropic two-phase medium results in two longitudinal waves and one shear wave, as predicted by Yakov Frenkel. The efficient and accurate numerical modeling of Biot’s equations of poroelasticity requires the knowledge of the exact stability conditions. This research presents the results of the von Neumann stability analysis of the discretized Biot’s equations and the discretized linear damped wave equation. The exact stability conditions for several implicit and explicit schemes are derived. Additionally, a multi-graphical processing units (GPU) numerical solver is developed to resolve the anisotropic elastodynamic Biot’s equations to simulate, in a few seconds, wave fields for spatial domains involving more than 4.5 billion grid cells. A comprehensive dimensional analysis is presented reducing the number of material parameters needed for the numerical experiments from ten to four. A dispersion analysis as a function of dimensionless parameters is performed leading to simple and transparent dispersion relations. The high efficiency of our numerical implementation makes it readily accessible to investigate three-dimensional and high-resolution scenarios of practical applications. As a part of the theory of non-linear instabilities, a new theory for earthquake nucleation is presented. The simplest visco-plastic or elasto-plastic rheology allows us to model spontaneous earthquake nucleation. By slowly increasing the stress in the medium, it reaches the yield surface, strain localization occurs resulting in the slow development of fractal shear bands. As time evolves, shear bands grow spontaneously, and stress drops take place in the medium. A stress drop corresponds to a particular new strain localization pattern which acts as seismic source and triggers seismic wave propagation. This new modeling approach is based on conservation laws without any experimentally derived constitutive relations

Redatuming and Quantifying Attenuation from Reflection Data Using the Marchenko Equation: A Novel Approach to Quantify Q-factor and Seismic Upscaling

Marchenko Imaging is a new technology in geophysics which enables to retrieve Green's functions at any point in the subsurface having only reflection data. This method is based on the extension of the 1D Gelfand-Levitan-Marchenko equation to a 3D medium. One of the assumptions of the Marchenko method is that the medium is lossless. If the lossy reflection response is used in the Marchenko scheme, some artefacts in the Green's functions as well as in the seismic image are present. One way to circumvent this assumption is to find a compensation parameter for the lossy reflection series so that the lossless Marchenko scheme can be applied. The main tasks of this thesis are to: [1] use the Marchenko equation to estimate the attenuation in the subsurface, [2] find a compensation parameter for the lossy reflection series so that the lossless Marchenko scheme can be applied, and [3] to create an upscaling method for wave propagation. The Artefact Removal Method was created which makes it possible to calculate an effective temporal Q-factor of the medium between a virtual source in the subsurface and receivers at the surface. This method is based on the minimization of the artefacts produced by the lossless Marchenko scheme. The minimization was performed in three ways: [1] in the space-time domain, [2] in the frequency domain and [3] to the scales of the wavelet transform applied to the artefacts. This method can also be used to find the layers with high attenuation. The upscaling method which can be used to construct macro-scale homogenized viscoelastic properties of the medium from the micro-scale properties of a heterogeneous medium was developed. This is done through linking the macro- and micro- scale Lippmann-Schwinger equations which describe the wave field and the strain field scattering in an inhomogeneous medium, respectively. In this thesis, the macro-scale homogenized viscoelastic properties were calculated by using the T-matrix Approach and the Generalized Dvorkin-Mavko Attenuation Model. All theoretical results are supported by synthetic 1D modeling. The theoretical part of the thesis and the general work flow can be used for a very complex medium